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. Abstract. General Schrodinger equation is considered with a central polynomial potential depend- 

' ing on 2q arbitrary coupling constants. Its exceptional solutions of the so called Magyari type (i.e., 

exact bound states proportional to a polynomial of degree -/V) are sought. In any spatial dimension 
D > 1, this problem leads to the Magyari's system of coupled polynomial constraints, and only 
purely numerical solutions seem available at a generic choice of q and N. Routinely, we solved the 
system by the construction of the Janet bases in a degree-reverse-lexicographical ordering, followed 
by their conversion into the pure lexicographical Grobner bases. For very large D we discovered 
that (a) the determination of the "acceptable" (which means, real) energies becomes extremely 
facilitated in this language; (b) the resulting univariate "secular" polynomial proved to factorize, 
■ utterly unexpectedly, in a fully non-numerical manner. This means that due to the use of the Janet 

\ bases we found a new exactly solvable class of models in quantum mechanics. 

1 Anharmonic Oscillators and the Problem of Their Solution 

Elementary Hamiltonian H = p 2 + q 2 + Xq i of the so called anharmonic oscillator in one spatial dimension 
D = 1 is an example which plays a key role in quantum theory and in many of its applications. We may 
recollect, for illustration, that small experimental irregularities in the vibrational spectra in atomic physics 
are currently being attributed to the quartic anharmonicity at a suitable and, if possible, reasonably small 
coupling constant, A = 0(1) [1]. For the fit of some experimental data of this type one may even employ 
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the two-parametric family of the Hamiltonians H = p 2 + q 2 + A q A + g q 6 [2], etc. In all these cases, 
sophisticated perturbation calculations are usually employed in order to achieve an agreement between 
experiment and theory (cf., again, ref. [1] and many other papers cited therein). 

In a mathematically more ambitious setting, Magyari [3] was probably the first who noticed that in one 
dimension, Schrodinger equation admits non-pcrturbativc, exceptional but exact bound-state solutions 

^(Magyari) ( g ) for 

any anharmonic potential of the following special polynomial form symmetric with 
respect to the origin, 

V^(r)^g r 2 +g 1 r i + ...+g 2q r^+ 2 , g 2q = 7 2 > , (1) 

provided only that its couplings gj satisfy certain q constraints. These constraints have the form of the 
system of coupled polynomial equations (their form will be displayed and discussed below) . Unfortunately, 
the achievement of the practical compatibility of the couplings with the Magyari's constraints requires 
the solution of his equations by a suitable more or less purely numerical technique. The corresponding 
algorithm is usually based on the use of Grobner bases [4] . The procedure is very standard and one would 
have no particular reason for its study in more detail in general. 

The first of the changes which proved relevant in this context appeared with the introduction of 
the higher-dimensional Schrodinger equations with polynomial interactions and with the Magyari-type 
solutions [5]. For all of these models, the Magyari-type equations become dependent on the dimension 
D > 1 playing the role of a new formal parameter. The freedom in its choice will prove most relevant 
in our present paper but in a historical perspective, it still took many years before this chance has been 
conceived and described in ref. [6] where the choice of the potential proved restricted, for purely technical 
reasons, to the "first nontrivial" polynomial of type (1) with exponent q = 2. 



The main inspiration of our present study of polynomial oscillators with q > 2 lies in the broad and not 
yet fully explored variety of the possibilities hidden in a consequent formal analysis of the Magyari-type 
equations. In this sense, the next decisive step has been made in refs. [7] where the Magyari's equations 
proved tractable in semi- and/or non- numerical manner at the first few lowest choices of the degree of the 
wave function, viz., at N = 1, N = 2, and N = 3. In these "trivial" cases, an overall tendency emerged of 
a distinct separation between the real (= physical) and complex (= apparently fully redundant) Magyari's 
couplings. This result offered an important hint for a more general analysis of the problem [8-10] and 
suggested the idea of using the Janet bases [11-13] in the similar cases. The simplicity of these low— N 
models enabled us to see that an overall and more systematic study should be directed towards the 
domain of the very large D ^> 1 (cf. also refs. [14-17] in this respect). 

In spite of the unique success of the mathematics of Magyari's nonlinear algebraic equations, a number 
of difficulties remained connected with their practical applications and applicability at the finite D. One 
of the key reasons (and differences from the harmonic oscillator and other exactly solvable models) is that 
the explicit construction of the Magyari's energies remains purely numerical. Indeed, these values (as well 
as the related couplings - we shall show some technical details below) must be computed as roots of a 
certain "secular" polynomial. This means that the difference between the variational, "generic N = oo" 
rule in Hilbert space seems only marginally simplified by the Magyari-type construction of any N 1 
bound state. 

The main purpose of refs. [6, 15] derived precisely from the latter point. Using the idea of perturbation 
expansions for Hamiltonians H = H^ q>N \D), these studies proceeded in two steps. Firstly, a zero-order 
approximations Hq (oo) have been constructed while, secondly, a series of corrections has been evalu- 
ated at each particular finite and fixed dimension D < oo. This opened the market for the constructions 
of the Hamiltonians Hq 9 ' N \oo) in systematic manner. 

In our present notation, the exact solvability of the zero-order Hamiltonians Hq Q ' N ^ emerged as an 
utterly unexpected result of our calculations at q = 1 in ref. [15], at q = 2 in ref. [6] and at q = 3 in ref. 
[17]. In what follows we intend to address the next, more sophisticated problems with q > 3. An emphasis 
is to be put on the vital role of the methods which were able to produce the necessary final results within 
the strict bound given by the not too fancy available computers. Hence, in what follows, the main emphasis 
will be laid upon the quality of the underlying software. Still, a more detailed introductory chapter is due 
first. 



2 The Derivation of the Magyari Equations 

2.1 Harmonic Oscillator with q — as a Methodical Guide 

The partial differential Schrodingcr equation for harmonic oscillator in D dimensions reads 

ti 2 1 \ 
-— A + -mO 2 \x\ 2 V(x) =e&(x) 
2m 2 / 



(2) 



and is solvable by the separation of variables in several systems of coordinates. The most common cartesian 
choice may be recommended for the first few lowest spatial dimensions D only [18]. In contrast, the 
separation in spherical system remains equally transparent at any D because it reduces eq. (2) to the 
same ordinary (so called radial) differential equation 
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ip(r) = Eil>(r) 



dr 2 

with r = \x\ e (0, oo), E = 2me/K 2 and t 
where L = 0, 1, . . .. At each L the energy levels are numbered by the second integer, 

E = E n x = u (2n + II + 3/2), n, L = 0, 1, . . . . 

The wave functions with quadratic A(r) = uir 2 /2 > and minimal N = n+ 1 in 
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ifi/h > 0. In this language we have I = II = L + (D — 3)/2 



(4) 



ip n ,L(r) = r 
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are proportional to an nth Laguerre polynomial [19]. In Hilbert space, their set is complete. 



2.2 q — 1 and Quasi-Exact (i.e., Incompletely Solvable) Sextic Oscillators 

An immediate partially or quasi-exactly solvable (QES) generalization of harmonic oscillators was dis- 
covered by Singh et al [2]. In this case one replaces 

oj — ► W{r) = a Q + a ir 2 , V (HO) (0) — ► G- X + G r 2 = U(r) (6) 

and gets the general sextic potential 

all three couplings of which arc simple functions of our initial three parameters and vice versa, 

g 2 =al>0, g 1 = 2a a 1 , g = 2al + G , G_ 1= 0. (8) 

The resulting Schrodingcr bound state problem cannot be solved in closed form. Nevertheless, we may 
postulate the polynomiality of the wave functions ip^ (r) for a finite multiplet (i.e., N— plet) of the 
wave functions. Under the specific constraint 

G = -ao-ai(4VV + 2^+l), N>1 (9) 
this AT— plet of polynomial solutions (5) is made exact by the choice of a WKB-like (i.e., quartic) exponent 

A(r) = \a r 2 + . (10) 

The ansatz (5) transforms then the differential Schrodingcr equation into a linear algebraic definition of 
the unknown iV— plet of coefficients h m . The solution is always obtained for a mere finite set of the levels 
n G (no,n\, . . . , njv— i). In contrast to the harmonic oscillator, the QES solvability is based on the L— 
and N— dependent constraint (9) so that, generically, the elementary QES multiplet exists in a single 
partial wave only. 



2.3 Magyari's QES Oscillators with q > 1 

The explicit energy formula (4) for harmonic oscillator was replaced by an implicit definition in the 
preceding paragraph which gives the sextic QES energies in the purely numerical form, viz., as zeros of 
the Singh's secular determinant of a certain tridiagonal N hy N matrix [20]. In this sense, Magyari [3] 
generalized the Singh's QES construction. In our present notation we may put, simply, 

V {q) {r) = U {q) {r) + r 2 [W (q) {r)] 2 , U (q) (r) = G r 2 + Gn A + ... + G q - X r 2q , 

W {q \r) = a a + a ir 2 + ... + a q r 2q (11) 

This formula re-parametrizes the polynomial (1) and specifies the one-to-one correspondence between the 
two sets of couplings, 

{9o, ■ ■ ■ ,92q} {G , ■ • . ,G 9 _i,a , ■ • ■ ,a q } 

where g 2q = a 2 , g 2q -i = g2q-i(a q , a q -i) = 2a,_i a q , . . . or, in opposite direction, a q = ^fg 2 l l = 7 > 0, 
otq-i = g2q-i/(2a q ) etc. 

At any q = 1, 2, . . ., equation (10) must be further modified, 

A (*>(r) = \a r 2 + + ... + ^-^a q r 2q+2 . (12) 

With a q > 0, one verifies that 

^{physical), r ^ _ e -A(«'(r)+0(l) i T ~> I 

which means that the correct bound-state ansatz 
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ip(r) = ]T fcW r 2n+i+1 exp -A^(r) (13) 



converts our radial equation (3) + (1) into an equivalent linear algebraic problem 
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with an asymmetric and non-square matrix 
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Its elements depend on the parameters in bilinear manner. 
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E - a (An + 21 + 3) 
-a 2 (4n + 21 - 1) + 2a ai - .gi, 



C„ = (2n + 2) (2n + 2£ + 3), 
i 1 ' = -ai (4n + 2^ + 1) + - 9o , 

. . . j 

A! 9 ' = -a g (4n + 2£ + 3 - 2q) + (a a g -i + ai« g -2 + • • ■ + a,-i a o) - ffg-i, 

71 = 0,1,... . 



(16) 



At any fixed and finite N = 1,2,... the non-square system (14) is an over-determined set of N + q linear 
equations for the N non- vanishing components of the vector h^ N \ At q = these equations degenerate 
back to the recurrences and define the harmonic oscillator states. At q = 1 we return to the sextic model 
where the "redundant" last row fixes one of the couplings and where we are left with a diagonalization of 
an N by N matrix which defines the N— plet of the real QES energies in principle. The situation is more 
complicated at q > 1. The counting of parameters and equations indicates that unless one broadens the 
class of potentials, only a very small multiplet of bound states may remain available in closed form [21]. 

Using an elementary change of variables, one may transform the decadic forces into their quartic 
equivalents etc. Paper [22] may be consulted for details which indicate that the study of any potential 
V(r) which is a polynomial in any rational power of the coordinate r may be replaced by the study of 
its present Magyari's or "canonical" QES representation V^ q '(r) at a suitable integer q. In addition, we 
shall also restrict our attention to the domain of large D. 

3 Magyari Equations at the Large Spatial Dimensions 

Up to now, our attention has been concentrated upon the structure of the QES wave functions. From the 
point of view of the evaluation of the energies, the main dividing line between the solvable and unsolvable 
spectra is in fact marked by the distinction between the closed q = formulae and their implicit QES 
form at q = 1. The transition to the next q = 2 may be perceived as merely technical. At all q > 1, 
the difficulties grow with N. In such a setting the emergence of certain simplifications at D 3> 1 may be 
crucial. 



3.1 An Exceptional, Decoupled Last-Row Constraint 

At any D, the last row in eq. (14) decouples from the rest of the system. At any q > 1 it may treated as 
a constraint which generalizes eq. (9), 

g q -i = -a q (An + 21 + 3 - 2q) + (a a q -i + u\u q -2 + ■ ■ ■ + a q -\ao) . (17) 

The insertion of this explicit definition of the coupling g q -\ simplifies the lowest diagonal in Q^ N \ 

A$ =47(^ + 9-71-!). (18) 



Since A^ +q _ 1 = we may drop the "hat" " and re-write cq. (14) in the more compact form where the 
size of the non-square matrix is merely (N + q — 1) by N, 

QlN] h (N) = Q _ (lg) 

This is the proper Magyari's system and it is merely solvable non-numerically in the simplest case with 
q = 0. No coupling is then fixed and the energies themselves are given by the explicit formula (17). Also 
the recurrences for coefficients of the wave functions may be solved in compact form. 

The next, q = 1 version of eq. (19) degenerates to the single, determinantal secular equation 

detQM=0. (20) 

Its solution is a purely numerical problem at all the larger N > 5. Of course, one coupling is fixed by eq. 
(17) and only the N— plet of energies must be calculated as represented by the real zeros of the single 
secular polynomial. 

At the larger exponents q > 2, some q mutually coupled N by N secular determinants must vanish 
simultaneously [14]. With an auxiliary abbreviation for the energy E = —g—i this means that at least 
one of the couplings is always energy-dependent and that its value must be determined numerically. In 
the other words, our non-square matrix = Q^(g-i, go, ■ ■ ■ ,g q -2) will annihilate the vector ft'"' if 
and only if all its q arguments arc determined in a deeply nonlinear and self-consistent, mostly purely 
numerical manner. 



3.2 Coupled Constraints at D ^> 1 

In our approach the guaranteed polynomiality of the wave functions will play a key role. One can say that 
in our original differential eq. (3) the numerical value of the spatial dimension D will be assumed large. 
No other simplifications will be assumed. 

In our problem with the old matrix elements 

C„ = (2n + 2) (2n + 2L + D), B n = -g_ x - a {An + 2L + D), 
A ( n ] = -g k -i -a k (4n + 2L + D- 2k) + {a a k _ x + ... + a fe _ia ) , (21) 
fc=l,2,...,g-l, n = 0,l,...,iV + g-2 

we shall preserve the dominant components of the matrix elements only, 

Cf = (2n + 2)D, 4°]=- fl _ 1 -a A = ~9k-x ~ «fc D , k < q 

(note that = is unchanged). Then we re-scale the coordinates and, hence, coefficients according 

to the rule 

h^= Pn /n n . (22) 

Simultaneously we have to replace the energies and couplings <7oj • ■ ■ j5?-2} by the new re-scaled 

parameters {s\, S2, . . . , s q } using the following linear recipe, 
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where we abbreviated 
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t = t(D) = (2 q+2 D q 1 



1/(9+1) 
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In the leading-order approximation this gives, finally, our Magyari equations in the compact form 
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which is to be studied in what follows. 
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(25) 



4 The Method of Solution of the D ^> 1 Magyari Equations 



4.1 Involutive Bases 

To solve polynomial systems (25) we shall construct for them the related Janet bases. Janet bases are 
typical representatives of general involutive bases of polynomial ideals [23] which are Grobnerian though, 
generally, redundant. However, just this redundancy of involutive bases makes the structural and combi- 
natorial information on polynomial and differential ideals and modules more accessible [24-27]. 

And as well as the reduced Grobner bases, the involutive bases can be used for solving polynomial 
systems with finitely many solutions that correspond to the zero-dimensional ideals [28]. For this purpose, 
a pure lexicographical monomial order seems best since it provides the completely triangular basis with 
sequentially eliminated variables starting from the highest one with respect to the order chosen [4]. 
However, computation of a lexicographical basis takes usually much more time than computation of a 
degree-reverse-lexicographical basis, first, and conversion of this basis into the lexicographical one, second. 

We use this two-step computational procedure in our study and solving (25). In doing so we shall deal 
with the minimal Janet bases [29] only. It is remarkable that a degree-reverse-lexicographical order is 
inherent in minimal Janet bases of zero-dimensional ideals. What follows from the demonstration in [26] 
is this inherence to Pommaret bases, and the fact proven in [30] that a minimal Janet basis is also a 
Pommarct basis whenever the latter exists, i.e., whenever it is finite. Zero-dimensional ideals always have 
finite Pommaret bases [31]. 

4.2 Janet Bases 

Below we use the definitions and notations from [29,30,32,33]: N is the set of non-negative integers; 
M = {xj 1 • • • xff | di G N} is the set of monomials in the polynomial ring R = K[xi, . . . , x n ] over zero 
characteristic field K; deg^w) is the degree of Xi in u G M; deg(-u) = Y^T=i ^ e Si( u ) ^ s the total degree of 
u; y is an admissible [4, 28] monomial ordering compatible with 

X\ >- X2 y ■ ■ ■ y X n . 

Divisibility of monomial v by monomial u will be written as u \ v. A divisor u of a monomial v is proper if 
deg(u) < deg(v). lm(/) and lt(f) denote, respectively, the leading monomial and the leading term of the 
polynomial / G R with respect to y. lm(F) denotes the leading monomial set for F, and Id(F) denotes 
the ideal in R generated by F. 

Let polynomial set F C R be finite and / G F. For each 1 < i < n we partition F into groups labeled 
by non- negative integers d\, . . . , g? ;: 

Mi, ■ ■ ■ , <U] = { f eF\d 3 = dc gj -(lm(/))> 1 < j < % }. 

A variable xi is (Janet) multiplicative for / G F if i = 1 and 

deg 1 (lm(/)) = maxjdeg^lm^)) \ g e F}, 

or if i > 1, / G [di, ■ ■ . , di-i) and 

degi(lm(/)) = maxjdeg^lm^)) | g € [di, . . . , di-i]}. 

If a variable is not multiplicative for / G F, it is nonmultiplicative for / and we write this as Xi G 
NMj(f, F). u G lm(i^) is a Janet divisor of w G M, if u \ w and monomial w/u contains only multiplicative 
variables for u. In this case we write u \j w. 

A finite polynomial set F is Janet autoreduced if each term in every / G F has no Janet divisors 
among lm(F) \ lm(/). A polynomial h G R is in the Janet normal form modulo F if every term in h has 
no J— divisors in 1m(F). We denote the Janet normal form of polynomial / modulo F by NFj(f, F). If 
the leading monomial lm(/) of / has no Janet divisors among elements in lm(F'), then we say that / is 
in the Janet head normal form modulo F and write / = HNFjlf, F). 

A Janet autoreduced set F is a Janet basis of Id(F) if any nonmultiplicative prolongation ( multipli- 
cation by a nonmultiplicative variable ) of any polynomial in F has vanishing Janet normal form modulo 
F: 

(V/ G F) (Vx G NMj(f, F)) [ NFj(f ■ x, F) = ] . (26) 

A Janet basis G of ideal Id(G) is minimal if for any other Janet basis F of the ideal the inclusion 
lm{G) C lm(_F) holds. A monic minimal Janet basis is uniquely defined by an ideal and a monomial 
order. In what follows we deal with the minimal Janet bases only and often omit the word "minimal" . 



4.3 Algorithm for Computing Janet Bases 



We present now the algorithm JanetBasis which is a special form of the general Gerdt-Blinkov algo- 
rithm [9, 29] for computing minimal involutive bases concretized for Janet division. This concretization 
in its more detailed form relied on the appropriate data structures - Janet trees - and is described in [32, 
33]. Note that, recently, the Gerdt-Blinkov algorithm in its form presented in [9,25] was implemented 
in Maple for both the polynomial and the linear differential ideals [12, 13]. 

To provide minimality of the output Janet basis [29] the intermediate data, i.e. initial polynomials 
and their prolongations and reductions, are partitioned into two subsets T and Q. Set T contains a part 
of the intermediate basis. Another part of the intermediate data contained in set Q also includes all the 
nonmultiplicative prolongations of polynomials in T which must be examined in accordance with the 
definition of Janet bases. 

To apply the involutive analogues of the Buchbcrgcr criteria and to avoid repeated prolongations we 
endow with every polynomial f E F the triple structure 



Here the ancestor of / is a polynomial g G F with u = fm(g) and such that u | lm(p). Moreover, 
if deg(u) < deg(lm(p)), then every variable occurring in the monomial lm(p)/u is nonmultiplicative 
for g. Besides, for the ancestor g the equality anc(g) = lm(g) must hold. These conditions mean that 
polynomial p was obtained from g, in the course of the below algorithm JanetBasis, by a sequence 
of nonmultiplicative prolongations. This tracking of the history in the algorithm allows one to use the 
involutive analogues of Buchberger's criteria to detect and avoid unnecessary reductions. 

The set vars contains those nonmultiplicative variables which have been already used in the algorithm 
for construction of nonmultiplicative prolongations. This set serves to prevent the repeated prolongations. 

After every insertion of a new element p in T all elements r G T such that lm(r) >- lm(p) arc moved 
from T to Q in line 13. Such a displacement provides minimality of the output basis. 

It should also be noted that for any triple p G T the set vars must always be a subset of the set of 
nonmultiplicative variables of pol(p). Line 21 controls this condition. 

The initialization step is done in lines 1-4. The subalgorithm JanetHeadReduce performs Janet 
reduction of the leading terms of polynomials in Q modulo polynomials in T. 

In the main loop 5-24 an element in Q is selected in line 6. The correctness of this selection strategy 
proved in [34]. In practice the cardinality Q at intermediate steps of the algorithm is rather large and 
easily runs up to hundreds and thousands. At the same time there may be different polynomials in Q with 
identical leading monomials. Therefore, the restriction in line 6 still admits some arbitrariness. In our 
implementation in [33] for the degree-reverse-lexicographical ordering a triple p G Q with the minimal 
dcg(lm(pol(p))) was chosen. In the case of several such polynomials in Q, the one with the minimal 
number of terms was picked up. 

Line 8 breaks computations in the case when inconsistency is revealed during the head term reduction 
in Q and returns the unit basis. In line 16 the tail Janet reduction is done, then the Janet reduced 
polynomial in p is inserted in T and all the higher ranked polynomials are moved to Q (loop 12-14). 
Actually this displacement takes place only if a polynomial in p has been subjected by the head term 
reduction in line 23. Otherwise, pol(p) >~ pol(r) holds for any r G T. The insertion of a new polynomial 
in T may generate new nonmultiplicative prolongations of elements in T which are added to Q in line 
20. To avoid repeated prolongations the set nmp(g) of Janet nonmultiplicative variables for q has been 
used to construct its prolongations is enlarged with x in line 21. 

The subalgorithm JanetHeadReduce computes the Janet head normal form of polynomials in Q 
modulo polynomials in T 

and invokes in line 6 subalgorithm HNFj (p,T) that does head reduction of a single polynomial p. 
For a head reducible input polynomial pol(/) the two involutive analogues of the Buchberger criteria [4] 
criteria are verified in line 8 of subalgorithm HNFj: 

— Criterion l(f,g) is true iff anc(/) • anc(g) | lm(pol(/)). 

— Criterion II(/, g) is true iff deg(lcm(anc(/) • anc(g))) < dcg(lm(pol(/)). 



P = {/, u, vars} 



such that 



pol(p) 
anc(p) 
nmp(p) 



/ is polynomial / itself, 

u is the leading monomial of a polynomial ancestor of / in F, 
vars is a (possible empty) subset of variables. 



Algorithm JanetBasis(F, -<) 

Input: FeR \ {0}, a finite polynomial set 

-<, an admissible ordering 
Output: G, a minimal Janet basis of Id(F) 

1: choose f € F with the lowest lm(/) w.r.t. y 

2: T := {/,lm(/),0} 

3: Q~{{q,\m(q),0}\qeF\{f}} 

4: Q :=JanetHeadReduce(<2,T) 

5: while Q / do 

6: choose p € Q such that lm(pol(p)) has no proper divisors among {lm(pol(q)) | q £ Q \ {p}} 
7: if lm(pol(p)) = 1 then 
8: return {1} 
9: else 

10: Q := Q \ {p} 

11: if lm(pol(p)) = anc(p) then 

12: for all {r £ T \ /m(pol(r)) >- lm(pol(p))} do 

13: Q:=QU{r}; T:=T\{r} 

14: od 

15: fi 

16: pol(p) := NFj(pol(p),T) 
17: fi 

18: T := T U {p} 

19: for all q G T and x G iVA/,/(pol(q), T) \ nmp(g) do 

20: Q — QUilpol^-s^nc^)^}} 

21: nmp(g) := nmp(g) n NMj(pol(q), T) U {x} 

22: od 

23: Q := JanetHeadReduce(Q, T) 
24: od 

25: return G := {pol(/) | / G T} 



Subalgorithm JanetHeadReduce(Q, T) 

Input: Q and T, sets of triples 

Output: Janet head reduced set Q modulo T 

1: S~Q 

2: Q:=0 

3: while S / do 

4: choose p £ S 

5: 5 := 5 \ {p} 

6: fc:=HNFj(p,T) 

7: if ft ^ then 

8: if lm(pol(p)) / lm(/i) then 

9: Q :=QU{h,lm(h),Q} 

10: else 
11: Q~QU{p} 
12: fi 
13: fi 
14: od 

15: return Q 



Subalgorithm HNFj(/,T) 
Input: / = {pol(/),anc(/),nmp(/)}, a triple 
T, a set of triples 

Output: ft = HNFj(po\(f),T), the Janet head normal form of the polynomial in / modulo polynomial set in T 
1: G := {pol( S ) | g G T} ' 

2: if lm(pol(/)) is involutively irreducible modulo G then 
3: return / 
4: else 

5: ft := pol(/) 

6: choose g G T such that lm(pol(<7)) ,/ lm(ft) 
7: if lm(ft) / anc(/) then 

8: if CriterionI(/, g) or CriterionII(/, g) then 

9: return 

10: fi 
11: else 

12: while ft 7^ and lm(ft) is L— reducible modulo G do 

13: choose q G G such that lm(g) |j lm(/i) 

14: /i := ft - 5 ■ It (ft) / It (g) 

15: od 

16: fi 

17: fi 

18: return ft 



If any of the two criteria is true, then HNF(pol(f),T) = [33]. Though as shown in [35] Criterion II 
does not fully replace the Buchbcrger chain criterion, in practice Criterion II works pretty well as our 
computer experiments demonstrate [33]. 

The last subalgorithm NFj performs the Janet tail reduction of a polynomial with irreducible leading 
term. It outputs the full Janet normal form NFj(f,T) of the input polynomial / modulo polynomial set 
containing in T. This subalgorithm is called in line 16 of the main algorithm JanetBasis and performs 
a chain of elementary involutive reductions until every term in the obtained polynomial becomes Janet 
irreducible modulo polynomials in T. 



NFj(/,T) 

Input: /, a polynomial such that / := HNFj(f,T); 

T, a set of triples 
Output: h — NFj(f,T), the full Janet normal form of ft 
modulo polynomial set in T 
1: G := {pol( S ) | g G T} 
2: h := / 

3: while h ^ and ft has a term t Janet reducible modulo G do 
4: choose g G G such that lm(g) \j t 
5: ft := ft -g ■ t/\t{g) 
6: od 

7: return ft 



It should be noted that both the full Janet normal form and the Janet head normal form are uniquely 
defined and, hence, uniquely computed by the above subalgorithms. This uniqueness is a consequence of a 
Janet divisor among the leading terms of polynomials in T at every step of intermediate computations [29] . 

4.4 Converting Bases and Finding Roots 

As we emphasize in the previous section, to find common roots of polynomials in a given system it is 
worthwhile to compute a pure lexicographical involutive or reduced Grobner basis. We do this computa- 
tion in the following three steps: 

1. Computation of a minimal degree- reverse- lexicographical Janet basis by the above described algo- 
rithm. 



2. Extraction from the Janet basis obtained the reduced Grobner basis. 

3. Conversion of the degree- reverse- lexicographical Grobner basis into the pure lexicographical one by 
the famous FGLM algorithm [36]. 

Step 2 is done immediately due to the history of prolongations stored in the polynomial triples (Sect. 
4.3). Since the reduced Grobner basis is a subset of the Janet basis computed [29] and this subset is 
irreducible with respect to the conventional (noninvolutive) reductions, a triple p = {/, u, vars} in the 
Janet basis contains an element / of the reduced Grobner basis if and only if lm(f) = u. This relation 
means that an element of the reduced Grobner basis is such an element in the Janet basis that it has 
no ancestors in the last basis. Indeed, the leading term of this element cannot be a prolongation of the 
leading term of other element in the basis. In addition to the use of criteria (Section 5.2) this is one more 
byproduct of the triple representation. 

The conversion of the degree-reverse-lexicographical Grobner basis extracted (old basis) into the pure 
lexicographical basis (new basis) is done as follows [36, 37]. First, a sequence of monomials is generated, 
starting from the least ones w.r.t. to the new ordering and then their normal forms are computed modulo 
the old basis until there appears a monomial whose normal form is a linear combination of normal forms of 
the preceding monomials. In this case, we add the polynomial given by this relation to the new basis. This 
process is continued by constructing other elements in the new basis by treatment of the next variables in 
accordance with the new monomial order. Computation of the normal form for a monomial is simplified 
if one takes into account the fact that the normal forms of all its proper divisors have been computed. 

A degree-reverse-lexicographical Grobner basis admits to find roots of the initial polynomial system 
by the sequential solving of univariate polynomial equations. Given a univariate polynomial, we tried first 
to factorize it and used the built-in factorization routines of computer algebra system Reduce 3.7 [38] for 
this purpose. If the factorization failed to give exact roots we used a special software package ROOTs 
written on the top of PARI-GP system [39] to find the roots numerically for the factors obtained. 



5 The Results for Polynomial Potentials with q < 3 
5.1 Sextic QES Oscillator with q — 1 and Any N 

Starting from the first nontrivial sextic-oscillator potential (7) with q = 1 and with the binding energies 
re-parametrized in accord with cq. (23) where si = s, 

E=\^LD + (64g 2 ) 1/4 VDs, 

2 y/92 

full attention must be paid to the selfconsistency problem represented by the set of equations (25). At 
every TV, its first nontrivial q = 1 version 



/ s 1 

A- 1 s 2 

N-2 s 3 



V 



2 s N — 1 
1 s 



( Po \ 

Pi 



PN-2 
\PN-lJ 



= 



(27) 



has the form of an asymmetric eigenvalue problem. In standard manner it leads to the secular equation 
(20) expressible as the following sequence of the polynomial conditions, 



4s = 0, 



10 s" 



9 = 0, 



20 s 3 + 64 s = 0, 



TV = 3, 

N = 4, 
N = 5, 



etc. By mathematical induction, all the infinite hierarchy of these equations has been recently derived 
and solved in ref. [15]. 

Quite remarkably, all of the real (i.e., "physical") energy roots s = s^' proved to be equal to integers. 
Moreover, all of them may be determined by the single and compact formula 



3 (i) 



-N- l + 2. ? , j = 1, 2, 



,N. 



(28) 



One imagines that all the coefficients p„ may be normalized to integers, 

= 1, N = 1, 

P^=P^=P^ = -P? = 1, N = 2, 

P? = ft W =pf =~pf ] =pf = pf = 1- P 1 1) = -pf ) =2, P f=0, N = 3, 

etc. 

The hrst result of our subsequent computations using the symbolic manipulation techniques proved 
equally encouraging since we succeeded in compactification of the set of the above recurrent solutions to 
the single leading-order form of the related wave functions, 

^)(r)=/ +1 (l + jY ' (i-j)' 'cxp^laor 2 -^ 4 ) , 

j = 1,2,. ..,7V. (29) 

A few more comments may be added. Firstly, the large and degenerate nodal zeros in eq. (29) are a mere 
artifact of the zero-order construction. This means that the apparently interesting exact summability of 
all the separate 0(r 2 //i) error terms is not too relevant, indeed. Although it leads to the zero-order nodes 
at r = 0(y/JI) — 0(D 1 ^ 4 ), these nodes have no real physical meaning. 

Secondly, the leading-order perturbative approximation provides a reliable information about the 
energies. They arc asymptotically degenerate, due to the large overall shift of the energy scale as explained 
in section 3.2. In addition, the next-order corrections may be easily obtained by the recipes of the textbook 
perturbation theory. As long as the coefficients p n are defined in integer arithmetics, the latter strategy 
gives, by construction, all the above-mentioned energy corrections without any rounding errors in a way 
outlined in more detail in ref. [15]. 

In the other words, we may say that formula (29) may either be truncated to its leading-order form 
exp (— \( 2 \r)) or, better, its full form may be used as a generating function which fa- 

cilitates the explicit evaluation of the coefficients p„ ' . In comparison, both the oversimplified harmonic 
oscillator and the q = 1 wave functions may be characterized by the similar coordinate dependence which 
becomes spurious (i.e., dependent on the selected normalization) everywhere beyond the perturbatively 
accessible domain of r. 

The energies specified by cq. (28) form an amazingly regular multiplct. A natural question arises 
whether a similar regularity could re-emerge at the larger integer indices q > 1. We are now going to 
demonstrate that in spite of the growth of the technical obstacles in dealing with the corresponding key 
equation (25), the answer is, definitely, affirmative. 



5.2 Decadic Oscillators with q — 2 and Any N 

The decadic anharmonic oscillator exhibits certain solvability features which motivated its deeper study 
in non-Hcrmitian context [40]. The changes of variables make this oscillator very closely related to the 
common quartic problem [14, 6]. Paying attention to the D ^> 1 domain and abbreviating the parameters 
si = s and S2 = t of the respective decadic-oscillator energy and coupling in eq. (23), we arrive at the 
four-diagonal version of our solvability condition (25) at q = 2, 

/ s 1 

t s 2 

N -1 t s 3 

7V-2 t s 4 



3 t s N -I 
2 t s 

V i t 

This is the first really nontrivial equation of the class (25) . In order to understand its algebraic structure 
in more detail, let us first choose the trivial case with N = 2 and imagine that the resulting problem 
(with pi due to the definition of N) may be solved by the determination of the unknown ratio of the 



Po 
Pi 



PN-2 
\PN-l 



(30) 



wave- function coefficients po/pi = —t from the last line, and by the subsequent elimination of t = 1/s 
using the first line. The insertion of these two quantities transforms the remaining middle line into the 
cubic algebraic equation s 3 = 1 with the single real root 8 = 1. 

The next equation at N = 3 is still worth mentioning because it shows that the strategy accepted in 
the previous step is not optimal. Indeed, the same elimination of P1/P2 = —t and of P0/P2 — (t 2 — s)/2 
from the third line leads to the apparently ugly result 

st 2 - s 2 -2t = 0, 
t 3 - 3st + 4 = 0. 

An alternative strategy starting from the elimination of po an( i P2 leads to the much more symmetric 
pair of the conditions 

t 2 - s 2 t + 2s = 0, 
s 2 - t 2 s + 2t = Q 

the respective pre- multiplication of which by t and s gives the difference < 3 = s 3 . This means that t = e s 
where the three eligible proportionality constants exist such that e = 1. Thus, our problem degenerates 
to a quadratic equation with the pair of the real roots s = t = s 1 - 1 ' 2 ^ such that 

s« = 2, s< 2 > = -1. (31) 

The "ugliness" of the procedure of elimination is inessential as long as we can produce the results by any 
"brute- force" symbolic manipulations on the computer [8]. 

Once we encountered the limitations of the naive algorithms, we were forced to pay attention to all 
the above-described sophistications of our algorithms. Fortunately this overall strategy proved successful. 
Using the methods described in preceding sections we revealed that the step-by-step elimination of the 
redundant unknowns gives the best form of the results when one uses the Janet bases. 

One of the main and most important byproducts of our approach is that the resulting final effective 
or "secular" polynomial equations for the single unknown quantity s depend in practice on its power s q+1 
only. The clear illustration is provided by the present q = 2 case at N = 3 giving the rule 

s 6 -7s 3 -8 = 0, N = 3. (32) 

This equation possesses the same complete set of the real roots (31) of course. Still, what is important 
is that any root r = Sq of eq. (32) itself still represents just the third (and in the more general cases, 
(q + l)st) power of the final relevant quantity with the physical meaning of the energy. Thus, we still 
have to solve the relation r = Sq where merely the value of r is known and where, therefore, at most one 
final parameter sq is real (= acceptable). 

One may conclude that the real energies of the "strongly spiked" decadic oscillator are very easily 
determined even without a detailed specification of an "optimal" elimination pattern, and that it is very 
easy to get rid of the redundant non-real roots sq at the very end of the algorithm. Of course, the numerous 
complex roots should not be discarded a priori as they might prove important in some other applications 
like a systematic computation of the corrections [15] which were not mentioned in our present paper at 
all. 

Our conclusions extracted at N = 3 are confirmed at the next N leading to the effective polynomial 
equation 

s 10 -27s 7 + 27s 4 - 729s = 0, TV = 4. 

Being tractable by our newly developed computer software and playing still the role of a test, it results 
in the set of the mere two real roots again, 

s (1) =3, s (2) =0, N = 4, q = 2. 

One finds that the q = 2 problem may be reduced to a single polynomial equation with ^ ^^j complex 

roots s at any N . The explicit calculations may be summarized in a statement that all the general physical 
(i.e., real) spectrum of energies proves to be quite rich and appears described by the closed and amazingly 
simple and transparent formula again, 



g(j) - 7v + 2 - 3j, j = 1,2,..., j m ax, jmax = entier 



N + 1 



(33) 



After one applies our Janet-basis procedure at the higher and higher dimensions N, one repeatedly arrives 
at the confirmation of the N— independent empirical observation (33) and extends it by another rule that 
at all the values of the dimension N, there exist only such real roots that = fw. This means that 
each "solvability admitting" real energy s requires, purely constructively, the choice of its own "solvability 
admitting" real coupling constant t. 



5.3 Oscillators with q — 3 and Their Solution at Any TV 

At q = 3 we have to solve the five-diagonal eq. (25), 
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which may be reduced, by means of the similar symbolic computations as above, to the single polynomial 
problem 

t 9 -I2t 5 - 64t = 

at N = 3, to the next similar condition 

t 16 - 68 1 12 - 442 < 8 - 50116 1 4 + 50625 = 
at N = 4, to the conditions of vanishing of the secular polynomial 

t 25 - 260 1 21 + 7280 1 17 - 1039040 1 13 - 152089600 1 9 + 2030239744 1 5 + 10485760000 1 
at N = 5, or to the perceivably longer equation 

t 36 - 777 1 32 + 1357 1 6 < 28 - 17 1 89 4 60 t 2i - 3513570690 < 20 - 
- 1 198527160446 1 16 + 103857100871252 1 12 + 873415814269404 1 8 + 
+74500845455535625 1 4 - 75476916312890625 = 

at N — 6 etc. These computations represent a difficult technical task but at the end they reveal again 
a clear pattern in the structure of the secular polynomials as well as in their solutions. One arrives at 
the similar final closed formulae as above. Now one only deals with more variables so that we need two 
indices to prescribe the complete classification scheme 
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t (j,k) = -/V-3 + 2j + 2fc, 
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(35) 



We may re-emphasize that all the real roots share the symmetry r = t but admit now a different second 
root s. The physical meaning of these roots is obvious. Thus, the energies of the oscillations in the 
polynomial well 

V(9=3,fc=l) ( r J = ar 2 + br 4 + ^ +gr U 

will be proportional to the roots r^' k > . After the change of variables, the roots will represent energies 
for the alternative, "charged" polynomial potentials 



etc [22]. 



6 The Results with q = 4 and q = 5 for N < iV, 



6.1 Non-Integer Roots Emerging at q — 4 and N < 6 



In our present formulation of the problem (25), we denote the descending diagonals as s m with m 
1, 2, 3, 4 and get the equation 
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(36) 



Its systematic solution does not parallel completely the above-described procedures. In fact, the reduction 
of the problem to the search for the roots of a single polynomial secular equation P(x) = (in the selected 
auxiliary variable x = — S4) enables us only to factorize P(x) on an extension of the domain of integers, 



P{x) = (x + 3) fox + l- y/fj (2x + 1 + v^) 
2 a; 2 — 3 x + 3 V5x + 18) (2 x 2 — 3 x — 3 VEx 



18 
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(-2V5 
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<V5x + 2x 2 ^ 



(-2 VE + 8 + 7 x - V5a; + 2i 2 j (2 Vf3 + 8 + 7 z + VEx + 2x 2 



2 x 2 + 2 x + 3 



V5)( 



2 + 2 x + 3 



\/5 + 3 - 3 x - \/5.t + 2 x 2 j f— V5 + 3 - 3 x + VEx + 2x 2 
(2y/E+8-3x + VEx + 2 x 2 ^j (-2 VE+8-3x- VEx + 2x 2 



From this lengthy formula it follows that we get 



3, 



(2) 



V5 + 1 



1.618, 



,(3) 



VE-i 



-0.618 



There only exist these three real roots S4 in this case. 

The similar computerized procedure gave us the real roots also at N = 5 and N = 6. The details may 
be found in ref. [17]. The inspection of these results leads to the conclusion that s 2 = S3 and Si = S4. 
We did not succeed in an application of our algorithms beyond N = 6 yet. The reason is that even the 
N = 5 version of eq. (36) in its reduction to the condition 



c 70 - 936 x 65 + 67116 a; 60 - 95924361 x 55 - 74979131949 : 



,50 



8568894879002 x 4 



... - 17459472274501870222336 x b + 142630535951654322176 = 

of the vanishing auxiliary polynomial required a fairly long computation for its (still closed and compact) 
symbolic-manipulation factorization summarized in Table 1 of ref. [17]. 



6.2 A Mind-Boggling Return of Integer Roots at q = 5 

N — 6 At q = 5 and TV = 6 the symbolic manipulations using the Grobner bases [4] generate the secular 
polynomial in x — S5 which has the slightly deterring form 

a 91 - 16120 a: 85 + 49490694 a; 79 - 286066906320 x 73 - 3553475147614293 x 67 - 

... - 319213100611990814833843025405983064064000000 a = . 

Fortunately, it proves proportional to the polynomial with the mere equidistant and simple real zeros, 

Pf'li) = x (a 2 - 1) (a 2 - 2 2 ) (a 2 - 3 2 ) (x 2 - 4 2 ) (a 2 - 5 2 ) . 

The rest of the secular polynomial is a product of the other two elementary and positive definite polyno- 
mial factors 

2 

p 2 ] 0) = II fc 2 ~ 3k x + 3fc2 ) ( x2 + 3fc2 ) + 3/c x + 3fc2 ) 

fe=i 

and 

5 

P 3 (6) = J| (a 2 - fca + /c 2 ) (a 2 + kx + k 2 ) , 

fe=i 

with another positive definite polynomial 

12 

p i 6) = n ( x2 - h x + cfe ) ( x2 + ^ * + c 
fe=i 

where the structure of the two series of coefficients (see their list in ref. [17]) is entirely enigmatic. 

The subsequent symbolic manipulations reveal a symmetry S2 = S4 and sj = S5 of all the real 
eigenvalues. In the N = 6 pattern summarized in ref. [17] we recognize a clear indication of a tendency 
of a return to the transparency of the q < 3 results which may be written and manipulated in integer 
arithmetics. For obtaining a deeper insight we must move to the higher N. 

N — 7 One should note that in spite of its utterly transparent form, the latter result required a fairly long 
computing time for its derivation. One encounters new technical challenges here. Indeed, the comparison 
of the N = 6 secular polynomial equation with its immediate N = 7 descendant 

x 127 - 60071 a 121 + 1021190617 a 115 - 11387407144495 a 109 - . . . + ca • 10 6 = 

shows that the last coefficient 

c= 125371220122726667620073789326658415654595883041274311330630729728 

fills now almost the whole line. This case failed to be tractable by our current computer code and offers 
the best illustration of the quick growth of the complexity of the q > 5 constructions with the growth of 
the QES dimension parameter N. 

Fortunately, we are still able to keep the trace of the pattern revealed at N = 6. Indeed, our new 
secular N — 7 polynomial factorizes again in the product of the four factors Pj(x), j — 1,2,3,4 where 
only the first one has the real zeros, 

P[ 7 \x) = p[ 6) (x)-(x 2 -Q 2 ) . 
The further three factors fit the structure of their respective predecessors very well, 
P 2 (7) (x) = P 2 (6) (a) • (a 2 - 9 a + 27) (x 2 + 27) (x 2 + 9 x + 27) 

and 

P 3 (7) = P 3 (6) ■ (a 2 - 6 a + 36) (a; 2 + 6 x + 36) 

while 

6 

P? ] = Pa ] ■ II ( x2 - h x + 9k) (x 2 + A a + g k ) . 

k=l 

The subscript-dependence of the new coefficients may be found in ref. [17] again. The key importance of 
the explicit knowledge of these coefficients lies in the possibility of a rigorous proof that the related roots 
are all complex and, hence, irrelevant from our present point of view. 



N — 8 and N — 9 The growth of the degree of our secular univariate polynomials makes it quite 
difficult to move too far with N. One may be more explicit in this respect: In place of the numerous 
irregularities encountered at q = 4. we may now be surprised by the re-emergence of the following closed 
and very transparent elementary formula for the q = 5 "energies" , 

s 5 €(-N+l,-N + 2,...,N-2,N-l). (37) 

which is valid again for any integer N in a way which parallels and complements the above-mentioned 
results which were available and published in our previous papers [15], [6] and [17] for the Magyari's 
D 1 potentials with q = 1, q = 2 and q — 3, respectively. In this context, their extension (37) is a 
brand new result which has not been published yet. Its unexplained equidistance property may be added 
to the list of the unresolved puzzles related to the Magyari equations. Indeed, the equidistance exemplified 
by eq. (37) would reflect a hidden algebra in linear cases but what is most intriguing here is the fact that 
the present Magyari equations are non-linear^. 

Another challenging feature of the problem lies in its exact solvability which is based on the factor- 
ization of polynomials of a very large degree V which grows, moreover, very quickly with N. Empirically, 
this degree may be even specified by the closed formula at q = 5 where T> = 3N 2 — 3N + 1 in a way 
illustrated by the next two explicit secular equations 

sf 9 - 186238 4 63 + 1176 88 1 3 1 99 4 57 ---- = 0, N = 8, (38) 

sf 7 - 502386 sf 1 + 94933635261 sf 5 - . . . = 0, N = 9. (39) 

On the basis of these observations we may conclude that an overall pattern of the smooth N— dependence 
of the equations survives, mutatis mutandis, smoothly the transition to the higher N. One can also prove 
(at least up to N < 9 at present) by construction that all the other factors of the secular polynomial have 
an elementary quadratic-polynomial form and remain positive for all the real "re-scaled energies" S5 . Their 
coefficients are elementary (we skip the examples here) so that the strict proof that they possess no real 
zeros is also elementary and very quick (one just evaluates the discriminants). A full parallelism between 
all N < 9 is achieved and we might conjecture, on this background, the possibility of its extension to all 
the non-negative integers TV. A strict proof of this conjecture could probably be based on mathematical 
induction but we do not feel it really urgent at the moment. 

7 Summary 

It is rather amusing to imagine that the majority of quantitative predictions in nuclear, atomic, molecular 
and condensed matter physics must rely on a more or less purely numerical model. The completely non- 
numerically tractable quantum systems are rare though, at the same time, useful and transparent (cf., 
e.g., the above-mentioned description of vibrations in molecules mimicked by harmonic oscillators). In 
our present paper we revealed that in the domain of the large spatial dimension D>1, the class of the 
exactly solvable models becomes, in a certain sense, broader. Thus, one might call all the polynomially 
anharmonic oscillators "asymptotically solvable" . 

This is an important and also not yet fully appreciated observation obtained due to the lasting 
advancement of the computer algebra and related software as described in more detail in Sections 4 and 
5. A fairly universal apparatus of these sections was reported in close connection with its application to 
our Magyari-type equations (25). 

In a certain perspective we found new closed solutions of Schrodinger equation with polynomial 
potentials in the domain of the large angular momenta £ 3> 1 where alternative techniques are also 
available (cf., e.g., their review [41] and/or very recent discussion [42]). Our results revealed the existence 
and provided the construction of certain fairly large multiplcts of "exceptional" £ ^> 1 bound states for a 
very broad class of polynomial oscillators. We believe that they might find an immediate application in 
some phenomcnological D 1 models. 

From the mathematical point of view, the most innovative and characteristic feature of our new 
D» 1 QES multiplets lies in the existence of the new closed and compact formulae for the QES energies 
and/or couplings at all N. For this reason, the corresponding partially solvable polynomial oscillator 

Hamiltonians H^ 1 '^ might even be understood as lying in the QES class as its new and fairly specific 
subclass. 

Due to an exceptional transparency of our constructions of HQ q ' N \ a facilitated return to the "more 
realistic" finite spatial dimensions D = 0(1) might prove tractable by perturbation techniques. Two 



reasons may be given in favor of such a strategy. First, due to the specific character of our present 
"unperturbed" spectra and eigenvectors, the perturbation algorithm might be implemented in integer 
arithmetics (i.e., without rounding errors) in a way outlined, preliminarily, in ref. [15] at q = 1. Second, the 
evaluation of the few lowest orders might suffice. This expectation follows from the enhanced flexibility of 
the available zero-order Hamiltonians. A priori, a better convergence of the corrections might be expected 
to result from a better quality of a "guaranteed smallncss" of the difference between a given Hamiltonian 
H at a finite D and one of its present D = oo QES approximants Hq. 
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